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Abstract 

A genetic algorithm approach suitable for solving multi-objective optimization problems is described and 
evaluated using a series ef aerodynamic shape optimization problems. Several new features includina two 
variations of a binning selection algorithm and a gene-space transformation procedure are included. The 
genetic algorithm is suitable for finding pareto optimal solutions in search spaces that are defined by any 
number of genes and that contain any number of local extrema. A new masking array capability is included 
allowing any gene or gene subset to be eliminated as decision variables from the design space. This 
allows determination of the effect of a single gene or gene subset on the pareto optimal solution. Results 
indicate that the genetic algorithm optimization approach is flexible in application and reliable. The binning 
selection algorithms generally provide pareto front quality enhancements and moderate convergence 
efficiency improvements for most of the problems solved. 
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set of scalar objective functions 
optimal set of F values 
scalar objective function 
maximum value of f 
rf h GA generation (active file) 
integer ranking array 

seed value used in the random number generator /7(0, 1) 

user-specified integer parameter controlling which selection algorithm is used 

user-specified integer parameter controlling the gene-space transformation option 

fixed number of chromosomes in each GA generation 

number of genes in each chromosome 

number of scalar objective functions 

user-specified vector with four elements that controls modification operator usage 
user-specified parameter controlling probability that a specific gene will be modified using 
perturbation mutation operator (0 < p < 1) 

user-specified parameter controlling probability that a specific gene will be modified using 
global mutation operator (0 < /£ < 1) 
basis vector matrix with elements n .- 

random number generator that returns a random value between 0 and 1 
unitary transformation matrix with elements i/. • 

/ A gene from the f h chromosome and the rf h GA generation 
f h chromosome from the GA generation 
user-specified maximum limit on the / h gene 
user-specified minimum limit on the f h gene 

user-specified parameter controlling the size of the perturbation mutations (0 < /3 < 1) 
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subscripts 

i 


J 

k 


gene or decision variable index 
chromosome index 
objective function index 


superscripts 

n GA generation or population index 

t temporary chromosome values obtained after selection but before modification 


Background 


Numerical methods for optimizing the performance of engineering problems have been studied for many 
years. Perhaps the most widely used general approach involves the computation of sensitivity gradients. 
These methods — called gradient methods — have been utilized to produce optimal engineering 
performance in a wide variety of different forms. The reliability and success of gradient methods generally 
requires a smooth design space and the existence of only a single global extremum or an initial guess 
close enough to the global extremum that will ensure proper convergence. 

In contrast to gradient-based methods, design space search methods such as genetic algorithms 
(GA) — also referred to as evolutionary algorithms (EA) — offer an alternative approach with several 
attractive features. The basic idea associated with the GA approach is to search for optimal solutions 
using an analogy to the theory of evolution. The problem to be optimized is parameterized into a set of 
decision variables or genes. Each set of genes that fully defines one design is called an individual or a 
chromosome. A set of chromosomes is called a population or a generation. Each complete design or 
chromosome is evaluated using a “biological-like” fitness function that determines survivability of that 
particular chromosome. For example, in aerospace applications, the genes may be a series of geometric 
parameters associated with an aerospace vehicle that is to be optimized for payload delivered to orbit, 
aerodynamic performance or structural weight. The fitness function takes as input all the geometric 
parameters and returns a numerical value for the fitness — the size of the payload, the aerodynamic 
performance or the structural weight. 

During solution advance (or “evolution” using GA terminology) each chromosome is ranked according to 
its fitness vector — one fitness value for each objective. The higher-ranking chromosomes are selected to 
continue to the next generation while the lower-ranking chromosomes are not selected at all. The newly 
selected chromosomes in the next generation are manipulated using various operators (combination, 
crossover or mutation) to create the final set of chromosomes for the new generation. These 
chromosomes are then evaluated for fitness and the process continues — iterating from generation to 
generation — until a suitable level of convergence is obtained or until a specified number of generations 
has been completed. 

Constraints can easily be included in the GA optimization approach either by direct inclusion into the 
objective function definition — a so-called penalty constraint — or, more commonly, by including one or 
more constraints into the fitness function evaluation procedure. For example, if a design violates a 
constraint, its fitness is set to zero, that is, it does not survive to the next generation. Because GA 
optimization is not a gradient-based optimization technique, it does not need sensitivity derivatives. It 
theoretically works well in non-smooth design spaces containing several or perhaps many local extrema. 

General GA details including descriptions of basic genetic algorithm concepts can be found in Goldberg , 1 
Davis , 2 and Beasley, et al . 34 Additional useful studies, which survey recent activities in the area of GA 
research including the presentation of model problems useful for evaluating GA performance are given in 
Deb , 5 Van Veldhuizen and Lamont 6 and Jimenez, et al . 7 
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A disadvantage of the GA approach is expense. In general, the number of function evaluations required 
for a GA optimization process to converge, exceeds the number of function evaluations required by a 
finite-difference-based gradient optimization (see the results presented in Obayashi and Tsukahara 8 and 
Bock 9 ). This situation is offset, to an extent, by the ease with which GAs can be implemented in parallel or 
distributed computing environments. 

Despite being relatively new, genetic algorithms have already been applied in many applications over a 
broad range of engineering fields. A brief survey of single discipline applications is presented in Holst 
and Pulliam , 10 and will not be discussed further in this report. 

Other applications involving GA search methods have been made in the area of multi-objective or multi- 
discipline optimization, that is, optimization problems in which two or more objectives are simultaneously 
and independently optimized. These methods, referred to as MOGA (multi-objective genetic algorithm) 
methods, are especially attractive because they offer the ability to directly compute so-called “pareto 
optimal sets” in a single computation instead of the limited single design point traditionally provided by 
other methods. 

The pareto optimal set or pareto front, as it is commonly called, includes optimal solutions for each of the 
individual objectives, as well as a range of tradeoff solutions in between, which are themselves optimal 
solutions. Providing a range of solutions to a multi-objective optimization problem is a powerful approach 
because it allows the designer to see the effect of decision variable variation on the design space in the 
form of optimal tradeoffs. Thus, the designer can choose individual objective weighting factors after their 
fun influence ts quantitatively known. 

In recognition of the importance of the MOGA approach, many theoretical developments have been 
published in recent years. In particular, Van Veldhuizen and Lamont 6 and Zitzler, et al . 11 present 
reasonably complete surveys of MOGA methodology with a special emphasis on how to compare GA 
performance from one algorithm variation to another. Other researchers, including Deb , 5 Deb, et al ., 12 
Van Veldhuizen , 13 Lohn, et al . 14 and Holst and Pulliam , 15 have developed and/or utilize a suite of test 
problems as a standard in evaluating MOGA convergence efficiency as well as the accuracy of the final 
pareto front. A key aspect of pareto front development is diversity. Does a particular MOGA produce 
good coverage over the entire pareto front or are some regions poorly resolved while other regions have 
high levels of undesirable clustering? This issue is studied in Laumanns, et al . 16 and Deb, et al . 17 

Still other MOGA issues are associated with archive strategies. As the pareto front develops, many 
solutions are found that lie on the pareto front that cannot be retained within the fixed population size of 
many schemes. How to retain or archive this information for the benefit of the MOGA while maintaining 
reasonable bounds on the archive file has been studied in Knowles and Crone . 18 - 19 One last example of 
MOGA research lies in the area of an interactive algorithm development. Parmee, et al . 20 present a MOGA 
approach that allows changes to be made in GA operators as well as in objective definitions as the MOGA 
advances from generation to generation. 

One area of MOGA research that is even more voluminous than the area of theoretical developments is 
that associated with GA multi-objective applications. Examples of multi-objective optimization applications 
include airfoil optimization by Marco, et al . 21 Naujoks, et al ., 22 Quagliarella and Della Cioppa , 23 Vicini and 
Quagliarella , 24 Hamalainen, et al . 25 and Epstein and Peigin , 26 missile aerodynamic shape optimization by 
Anderson, et al ., 27 and wing optimization by Anderson and Gebert , 28 Sasaki, et al ., 29 Oyama , 30 Obayashi, 
et al ., 31 and Ng, et al . 32 

Additional MOGA examples in the area of turbomachinery optimization include rocket engine turbopump 
design by Oyama and Liou 33 and compressor blade design by Benni 34 and Oyama and Liou . 35 In some of 
these examples the multiple objectives were obtained by considering two different aerodynamic design 


3 



points. In others the multiple objectives involved different disciplines including aerodynamics, structures, 
controls and/or electromagnetics. 

One last area of GA research that bears mention is In the area of hybrid methods, the utilization of GA 
optimization in conjunction with another type of optimizer. This is an attractive area of research because 
GA methods are particularly good at finding a global extrema, but not particularly good in converging the 
optimal solution to tight levels of accuracy. Briefly, several examples where hybrid approaches have been 
utilized include Rai 36 and Madavan 37 for single objective problems and Giotis, et al., 38 Tursi, 39 Vicini and 
Quagliarelia 40 and Brown and Smith 41 for multi-objective problems. 


Various definitions and the multi-objective genetic algorithm used in the present study are described 
next. Details associated with each of the operators, including selection, passthrough, random average 
crossover, perturbation mutation and mutation are presented. MOGA effectiveness is then evaluated 
using several three-dimensional aerodynamic shape optimization problems. 


Problem Statement: Single Objective Optimization 

A single-objective optimization problem can be stated as follows: Let / be a scalar function of N G 
independent variables, x n defined on some domain Q 

f= f{X) = f{x v --,x h -,x Ng ) (la) 

In this notation X is the vector of design space decision variables. The maximum value of f , indicated by 
/*, is obtained by finding the values of X = X r such that 1, 

f = max{/} = f(X*) = f{x;,-,x';,-,x„ 0 ) (1 b) 

The above maximization operation is subject to N co conditions or constraints indicated by 

c„(X)<0 n = \2,-,A/ co (1c) 

The constraints placed on the decision variable vector X by Eqs. (1c) essentially serve to limit the design 
space within Q for which the optimal solution is sought. 


Problem Statement: Multi-Objective Optimization 

For optimization problems involving more than one objective, which are simultaneously and 
independently optimized, the situation is more difficult. This is because each objective must play a role in 
determining the optimal solution. In the optimization process, conflicts might arise among the various 
objective functions, that is, the optima! values of each individual objective, in general, will not occur for the 
same decision variable vector, X . As a result, the “optimal solution” for a multi-objective optimization 
problem is typically a range or a set of solutions, which represent tradeoffs in objective space. 

To determine when one solution is better than another for multi-objective optimization problems the 
concept of dominance is utilized. 1 A vector U = U(^, •••,£//, -•-,%) is said to dominate another vector 
V = V(ir 1 , -,v N ) if and only if U/ > v, for all / and there exists at least one value of / such that u, > v r A 
vector defined on some domain Q that is not dominated by any other vector defined on O is said to be 

non-dominated on Q. 


r For the purpose of simplifying the discussion of algorithmic details, maximization is generally assumed. 
The logic for minimization is a straightforward modification and will not be discussed. 



A multi-objective optimization problem can be stated as follows: Let F be a set of N 0 scalar objective 
functions, f k , each dependent upon the same decision variable vector, X, which is defined on some 
design space Q 

F = F (2a) 

As above, the decision variable vector X consists of N e independent components. The multi-objective 
optimization problem involves finding the set of X = X* that produces non-dominated values for F = F* on 
£2. This set of values F* is called the pareto optimal set or the pareto front. 

For each F the constraints 

^(X)<0 n=\2,--,N co (2b) 

must be satisfied. Existence of these constraints serves to limit or reduce the size of Q, for which the 
optimal solution is sought. 


Genetic Algorithm 

The genetic algorithm optimization procedure utilized to solve the multi-objective optimization problem, as 
described by Eqs. (2), is now presented, it is closely related to the MOGA optimization procedure 
presented in Holst and Pulliam. 15 As mentioned in the introduction, the general idea behind GA 
optimization is to discretely describe the design space using a number of decision variables, x,. In GA 
parlance these parameters are called genes, and the / subscript refers to the gene number. Each set of 
genes that leads to the complete specification of an individual design, that is, each decision variable 
vector, X , is called a chromosome and is indicated by 


X'? 

A / 





(3) 


where the / subscript, added to X , identifies the chromosome number. In addition, the / subscript has 
been added to each gene, to indicate its chromosome membership. The n superscript has been added 
to indicate the GA generation number, which is iteratively advanced as the solution converges. Thus.X" is 

the f' chromosome for the /7 th generation that consists of N a genes. 

For aerodynamic shape optimization problems, the design space genes are typically a series of geometric 
parameters, for example, airfoil thickness and camber and/or wing sweep, twist and taper. For many GA 
applications genes are computationally represented using bit strings and the operators used to 
manipulate them are designed to accommodate bit string data. In the present approach, following the 
arguments of Oyama, 42 Houck, et al. 43 and Michalewicz, 44 real-number encoding is used to represent all 
genes. 

The key reason for using real number encoding is that it has been shown to be more efficient in terms of 
CPU time relative to binary encoded GAs. 44 In addition, real numbers are used for all genes in the present 
implementation because the present engineering application involves decision variables that are best 
described using real numbers, for example, the geometric parameters in aerodynamic shape optimization. 
Thus, using real number encoding eliminates the need for binary-to-real number conversions. 
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Initialization 


unce tne Design space has boon defined in terms of a sat of real-number genes, tha naxt stap is to form 
an initial generation, G°, represented by 


G°=(X?,--,X°v,X^) (4) 

where N c is the total number of chromosomes. Each gene within each chromosome is assigned an initial 
numerical value using a process that randomly chooses numbers between fixed user-specified limits. For 
example, the / <h gene in an arbitrary chromosome is initialized using 

Xj — /F\0, 1)(-^ maX/ — -*rnin, ) ''■min, (®) 

where x m3X . and x min/ are the upper and lower limits for the / h gene, respectively, and /7(0, 1) is a random 
number generator that delivers an arbitrary numerical value between 0.0 and 1 .0. 

The random number generator used in the present study requires an integer input — a seed value. If the 
integer is positive, the next number in the current random number sequence is returned. If the integer is 
negative, the random number sequence is reset to a new sequence. Utilization of the same negative 
seed value will always reset the random number generator to the same random sequence. Each new 
solution begins by resetting the random number generator using a single call to /?(0,1) with a negative 
seed value — ISEED. All other calls to /j\0,1) during that solution use a positive seed value. Thus, a 
solution can be repeated by simply using the same initial seed value or rerun to determine statistical 
variation by using a different initial seed value. 

Fitness evaluation 


After a generation is established — either the initial generation or any succeeding generation in the 
evolutionary process — fitness values, Ff, are computed for each chromosome using a suitable function 

evaluation. This is analogous to the objective function evaluation in gradient methods and is represented 
using 


F/=F(X?) (6) 

For example, fora multi-discipline optimization (MDO) problem involving the simultaneous maximization of 
two separate and distinct objectives, 4 anc * 4- ^e fitness evaluation represented by Eq. (6) consists of 
the following 

f = {(*/) 

? = 4(X/) 

where, for example, the first function jf might be the aerodynamic drag of an aerospace vehicle 
(constrained to fly at fixed lift) and the second function 4 might be the structural mass of the same vehicle. 
Once all the genes in a specific chromosome have been specified, j( can be evaluated using a suitable 
CFD solver to obtain the drag and £ can be evaluated using a suitable structural analysis routine to obtain 
the structural mass. In this case, of course, the optimization problem would be one of minimization. 
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Ranking 


The purpose of the ranking operation is to determine a set of integer values for the ranking array, !R. One 
integer value is required for each chromosome. The ranking algorithm is quite different for multi-objective 
optimization problems relative to single-objective problems. As such, both situations will be discussed in 
this section. The ranking array values — once determined — are then used in the GA selection process. 

Single Objective Optimization — The GA ranking algorithm for single objective optimization problems is 
quite simple. Whichever chromosome has the highest fitness is ranked number one (//?=1), whichever 
has the second highest fitness is ranked number two (//?= 2), and so on. This ranking algorithm, more 
formally stated, is given by 


ic-\ 

<fjj) ic- ic+~\ jj = \ - 
!R” = ic 


N r 


j- \ -,N c 


(7) 


where j and jj are special counters that range over all N c chromosomes in the current population or 
generation level. 

Multi-Objective Optimization — For multi-objective optimization problems the ranking procedure is more 
complex and utilizes the concept of dominance, which was defined previously. 

A chromosome with a fitness vector F that is not dominated by any other fitness vector within the design 
space is said to be a non-dominated chromosome. The optimal solution set F* or pareto front includes all 
solutions with fitness vectors that are non-dominated and that satisfy the constraints given by Eqs. (2b). 
The numerical approximation to the pareto front must involve a suitably complete set of discrete solutions 
so as to describe the. optimal values of each individual objective — these are the pareto front 
endpoints — as well as the non-dominated tradeoff solutions in between the optimal values associated 
with each of the individual objectives. 

With the concept of dominance in hand the actual ranking process can now be presented. There are a 
number of ranking procedures available for use in multi-objective optimization. The one used in the 
present study is called Goldberg ranking.' After all of the fitness values have been determined, each 
chromosome is checked for dominance. Those chromosomes that are non-dominated are given a number 
one ranking (/Ft= 1) and then temporarily removed from the current generation. The remaining 
chromosomes that are non-dominated are given a number two ranking (//?=2) and then temporarily 
removed. This process continues until each chromosome has an integer value for the ranking array, //?. In 
general, with this approach, the number of different integer values contained within the ranking array will 
be small, at least small in comparison to N c , as many chromosomes will be ranked near the top with a value 
of 1 , 2 or 3. 


The above procedure, by itself, represents a legitimate algorithm for the ranking operation. However, 
there is a refinement that provides a considerable enhancement to MOGA convergence. Before 
describing this refinement, two additional concepts must be defined — the active file and the accumulation 
file. 


The active file is simply a specific name used for the current generation of chromosomes, that is 


g" = (x(V^x;,-,x£j 
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is the /7 th generation active file for the GA iteration process. As the GA solution evolves, the active file is 
always fixed in size at N c chromosomes. The first N a elements in the active file — for the present 
implementation — always contain the chromosomes that posses the maximum fitness values for each of 
the N 0 objectives. 


The accumulation file is the list of all non-dominated chromosomes that have been found from all 
generations combined during the current MOGA iteration. As the MOGA solution evolves, the 
accumulation file typically grows in size with more chromosomes being added after each new generation. 
As new non-dominated chromosomes are added to the accumulation file, old chromosomes that lose their 
dominance must be deleted, thus ensuring that the accumulation file contains nothing but number-one 
ranked chromosomes. Because each chromosome stored in the /7 th generation accumulation file is non- 


it s o rvs s to define the psrsto front or 


at least, the rj generation approximation to the pareto 
front. The use of an accumulation file makes sense onlv when Nr > 2. 


The multi-objective ranking procedure described above utilizes only the chromosomes in the active file, 
that is, each chromosome in the active file is ranked relative to the other chromosomes in the active file. 
But this philosophy potentially wastes a plethora of information because not every number-one ranked 
chromosome for multi-objective optimization problems can be retained from one generation to the next in 
the active file. 

This is where the refinement in the ranking routine enters. Once each chromosome in the active file is 
ranked using the standard routine, an additional test is performed to see if any number-one ranked 
chromosomes in the active file are dominated by any of the chromosomes in the accumulation file. If this is 
the case, then the ranking number associated with the newly ranked chromosome in the active file is 
decremented by one. This ensures that the ranking routine produces number one rankings that are 
number one in a global sense. 

Selection 


After the ranking array is established, with or without the accumulation file option, the GA algorithm passes 
to the selection process to determine which chromosomes will continue on to the (/7+1) st generation and 
which will not. In the present algorithm implementation, four different selection operators have been 
coded: (1) a simple technique called greedy selection, (2) a traditional tournament selection algorithm, (3) 
a bin selection algorithm based on the computation of arc-length along the pareto front and (4) a second 
bin selection algorithm based on subdividing the design space into a matrix of boxes. The first two 
selection algorithms can be used for one or more objectives; the fourth can be used for two or more 
objectives, and the third inherently requires two objectives. The selection algorithm actually used is 
controlled by a user-input parameter called ISLECT, which will be defined shortly. 

Regardless of the selection algorithm chosen, the selection process picks established 
chromosomes — either from the /7 th generation active file, G", or from the accumulation file — and then 
stores the results in a temporary holding array G'. The chromosomes stored in G' are then used by the 
subsequent modification operators to create G n+1 , which will be discussed shortly. For the present study, 
results are computed using only the first, third and fourth selection algorithms, and thus, only these will be 
described herein. 
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Greedy selection— This selection algorithm is quite simple, it selects chromosomes from the / 7 th 
generation active file, that is, from Xy. it is implemented using the following 


rrr= J \ 


if (/Ay 7 < if] then 
yj - y n . 

/77 = / 77+1 


-f=\N c 


it=\N c 


if(/77> A^istop 
endif 


( 9 ) 


where each selected chromosome X f m is placed in a temporary holding array indicated by 

Note how the highest ranked individuals in the /7 th generation are selected multiple times, individuals with 
average ranking are selected a small number of times and individuals with the lowest rankings are not 
selected at all. This biasing toward individuals with the highest rankings is a key element in any GA. This 
baseline variation of greedy selection is implemented by setting ISELECT = 1. 

Tournament selection — The tournament selection algorithm is quite simple and is used widely — one 
variation or another — in many GA optimization applications. The present variation is implemented as 
follows. Using a random process, three chromosomes are chosen from the / 7 th generation active file, that 
is, from X'J. The ranking array values, //?, are then compared. The chromosome with the highest ranked 

array value is placed into a temporary holding array, G'\ In case of ties the chromosome selected first is 
chosen. Next, three new chromosomes are chosen from Xj , and their ranking array values are compared. 

As before, the chromosome with the highest ranking is added to G'. This process continues until N c 
chromosomes have been selected from Xj. At the conclusion of the tournament selection process, 
some of the original chromosomes within Xj may not have been selected even once, while other 

chromosomes may have been selected multiple times. This selection algorithm is implemented by setting 
ISELECT = 2. 

Bin selection (version 1) — The first bin selection algorithm is different from greedy and tournament 
selection, as implemented here, in two general ways. First, the bin selection algorithm chooses its 
chromosomes from the accumulation file, not the active file. Second, the bin selection algorithm divides 
the distance along the pareto front into N bin equal segments or “bins” using an arc-length computation 
and then selects an equal number of chromosomes from each segment. This ensures an equal 
distribution of selected points along the pareto front. The parameter N bin is a user-controlled parameter 
typically equal to 5 or 1 0. 

More precisely, this bin selection algorithm (version 1 ) is implemented as follows: 

(1) initiate optimization using greedy selection and proceed until the iotai number of points on the pareto 
front equals or exceeds N tot , which is a user controlled parameter that is typically several times the value of 

^birr 


(2) The arc-length along the pareto front is computed then divided into N bjr , equal arc-length segments or 
bins. Each existing pareto front solution point is assigned to a bin according to its arc-length value. 
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Depending upon the distribution of points along the pareto front, some bins may be heaviiy populated 
and others lightly populated. 

(3) The selection process randomly chooses a chromosome from the accumulation file keeping track of 
which bin it was chosen from. When the number of chromosomes chosen from a particular bin equals N avg 

defined by 


^avg ^ 'bin 

further selections from this bin are blocked. The process continues until N c chromosomes have been 
selected. For situations when N c is not evenly divisible by N bin , the value of N avg is incremented by one. 

Thus, certain bins may supply one more element to the next generation than other bins, out generally, the 
distribution wiii be reasonably weii distributed across the entire pareto front. This baseline variation of the 
bin selection algorithm (version 1) is implemented by setting ISELECT = -3. 

This selection algorithm does not guarantee that the pareto front endpoints will be selected and carried 
forward to the modification process or passed through to the next generation. Since endpoint retention 
can be an attractive feature for selection, a variation of this selection algorithm that first selects the pareto 
front endpoints and then proceeds with the standard selection algorithm is also available and is 
implemented when ISELECT = 3. 

Since the present bin selection algorithm requires the computation of arc-length along the pareto front, it 
inherently requires optimization problems with two objectives, that is, the pareto fronts must be a curve in 
two-dimensional space. Another bin selection option, which works for two or more objectives, is available, 
and will be discussed next. 

Bin selection (version 21 — The second bin selection algorithm is similar to the first in philosophy and has 
virtually identical convergence properties, but differs in how the bins are constructed. It does not use an 
arc-length computation along the pareto front. Instead, the design space is divided into a series of equally 
sized boxes or bins, each with /^-dimensions. This approach is applicable for any number of objectives 
(greater than one) and thus is more general than the first bin selection algorithm. 

More precisely, this bin selection algorithm (version 2) is implemented as follows: 

(1) Initiate the optimization using greedy selection and proceed until the total number of points on the 
pareto front equals or exceeds N tot — defined the same as in bin selection algorithm version 1 . 

(2) The design space bins are constructed next by subdividing each objective dimension in the design 
space into equal segments — where is a user-controlled parameter. Thus, the design space is 

divided into a total of (M seg ) No bins. For example, when U seg = 5 and = 2 , the design space is divided 
into 25 bins, each with the same rectangular shape. 

(3) Each existing element on the pareto front, that is, each element in the accumulation file, is categorized 
as to which bin it occupies. Many of the design space bins will be empty, because they lie beyond the 
pareto front or because they occupy an interior portion of the design space. Some of the bins that contain 
pareto front elements may have many entries — some may have few. This, of course, depends on how the 
pareto front is populated and how the pareto front slices through the matrix of bins that have been 
created. The number of bins that have at least a single entry is counted. That resulting value is then stored 
in N bin . 
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(4) The selection process randomly chooses a chromosome from the accumulation file keeping track of 
which bin it was chosen from. When the number of chromosomes chosen from a particular bin equals 
ined the same as above — further selections from this bin are blocked. The process continues 


avg 

until N c chromosomes have been selected. As above, for situations when N c is not evenly divisible by 
N bjn , the value of N 3vg is incremented by one. Thus, certain bins may supply one more element to the 

next generation than other bins, but generally, the distribution will be reasonably well distributed across 
the entire pareto front. This baseline variation of the bin selection algorithm (version 2) is implemented by 
setting ISELECT = -4 . 


As with bin selection algorithm (version 1), this selection algorithm does not guarantee that the pareto 
front endpoints will be selected and carried forward to the modification process or passed through to the 
next generation. A variation of this algorithm that first seiecis the pareto front endpoints and then 
proceeds with the baseline selection algorithm is also available and is implemented when ISELECT = 4. 

Modification Operators 


P Vector — After the new chromosomes have been selected and placed in the temporary holding array, 
they must be modified using one of several modification operators to obtain the (/7 + 1) st generation of 

chromosomes, G n+1 . In the present implementation four modification operators are used — passthrough, 
random average crossover, perturbation mutation and mutation. How many chromosomes are modified 
with each operator is controlled by the ^vector, which consists of four parameters — p B , p A , p P , p M . 
Each element of the P vector controls one modification operator. The value of each P vector element 
ranges from 0 to 1 .0, and, for consistency, the sum of all four elements must always equal one. A P vector 
equal to 0.1 , 0.3, 0.3, 0.3, for example, will cause the first 10 percent of the chromosomes to be modified 
using the passthrough operator, the next 30 percent to be modified using random average crossover, the 
next 30 percent to be modified using perturbation mutation and the last 30 percent to be modified using 
mutation. That is, p B = 0.1, p A = 0.3, p P =0.3, and p M = Q. 3. 


The passthrough operator is always performed first. After passthrough is complete, the implementation 
order of the remaining operators is immaterial. Once all values of G n ^ have been established, the 
algorithm proceeds to fitness evaluation, ranking and onto succeeding generations until the optimization 
is sufficiently converged. 

Passthrough — The simplest modification operator used in the present GA is “passthrough.” As the name 
implies, a certain number of chromosomes with the highest individual fitness values are simply “passed 
through” to the next generation from G / to G ^ 1 without modification. The number of chromosomes that 
are passed through to the next generation is controlled by the first parameter in the P vector, p B . The 

passthrough operator is always performed on the first p B N c chromosomes in G / . Care must be taken 
when choosing p B and N c so that p B N c >N 0 . If this is done and if a selection option which retains end 
points is used, the chromosomes with the highest individual fitness values will always be passed through 
to the next generation, thus guaranteeing that none of the individual maximum fitness values will ever 
decrease during GA iteration. 

Random average crossover — The next GA modification operator is called random average crossover and 
is implemented by first selecting two random chromosomes Xy, and Xy 2 from G / . Next, the two selected 
chromosomes are combined on a gene-by-gene basis using the following formula: 


+ 4 / 2 ) '=i. 2,-, 


N n 
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where xj’J' 1 corresponds to the / <h gene in the / b chromosome associated with G^ 1 and xf y1 and xf y2 

correspond to the / <h genes from the randomly chosen chromosomes X y1 and X y2 . The number of 

chromosomes modified using the random average crossover operator is determined by the parameter 
p A — the second element in the P vector. 

Perturbation mutation — The next GA modification operator is called perturbation mutation and is 
implemented by first selecting a random chromosome X y - from G'\ Next, a probability test is performed for 

each gene in the selected chromosome involving a call to the random number generator /7(0,1) . If the 
returned random number is less than p the gene is modified using 

< +1 = ^ + (x maX; - x min< )[AK0,1) - 0-5]/3 (11) 

where P is a user-specified tolerance that controls the size of the perturbation mutation, and p is a user- 
specified control parameter that statistically controls the number of genes that are modified. For sensible 
results the values of p and p must be between 0 and 1 .0. 

Because this operator can cause the value of a particular gene to exceed one of its constraints ( x maX/ or 
x mjn ), checks are required to prevent this. The number of chromosomes modified using the perturbation 
mutation operator is determined by the parameter p P — the third element in the P vector. 

Mutation — The last GA modification operator used in the present study is called mutation and is 
implemented similarly to the perturbation mutation operator. First, a random chromosome X y - is chosen 

from G / . Next, a probability test is performed for each gene x£ y in the selected chromosome involving a 
cal! to the random number generator ^0,1). if the returned random number is less than p% the gene is 
given a completely different value using 

~ (''max, ~ Tnin, "0 " r ''min, 02) 

The parameter p^ is a user-specified control parameter that statistically controls the number of genes that 
are modified. For sensible results p~ must be between 0.0 and 1.0. The number of chromosomes 
modified using the mutation operator is determined by the parameter p M — the fourth element in the P 
vector. 

Gene-space transformation — An option for accelerating GA convergence for multi-objective optimization 
problems involving a gene-space transformation is described next. This option, introduced in Ref. 15, 
worked amazingly well for model problems with simple non-convoluted pareto fronts, but performed 
poorly for model problems in which the pareto front was convoluted. Thus, it is of interest to see how the 
gene-space transformation procedure will perform on a real-world problem in which the shape of the 
pareto front — in gene space — is unknown. 

The transformation algorithm is outlined as follows: 

(1) Transform G' using the gene-space transformation procedure. 

(2) Perform modifications on the transformed chromosomes according to the user-input P vector 
values. 

(3) Using the inverse of the original transformation the modified chromosomes are transformed back 
to obtain G"" 1 . 
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This gene-space transformation procedure, which can be viewed as a gene-space rotation of coordinates, 
causes a linear coupling between each of the genes, and thus, affects how they are changed in the 
modification process. In some cases the gene-space transformation procedure can significantly improve 
GA convergence. The transformation procedure, which theoretically works for any number of objectives, 
will now be presented for the two-objective case, that is, for N a -2. 

The idea behind the transformation procedure is to perform a rotation of coordinates in gene space using 
an angle-preserving, length-preserving orthogonal transformation. For this purpose a simplified version of 
the Gram-Schmidt orthogonalization is used. 45 The current set of /7 th generation chromosomes (those that 
have been newly selected and placed in the temporary holding array GO can be written as 

X\M- \ 

■ w I 


Hg.Hc) 
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G' = 


1.2 


* 2,1 


*Afe,1 


where each column is one of the chromosomes in the active file holding array. It is this matrix that is to be 
transformed using 


UG / = (13) 

where U is a unitary matrix of rank N G that needs to be constructed and G' 7 is the resulting transformed 
matrix. To construct U a set of basis vectors R needs to be established. The elements of R are defined by 

tjj - X 2 - X., for /=1 (1 4a) 

and 

fl if /=/ 

Ho » /*/ ,or '- 2 <14b) 

The simple choices made for r, : when /> 2 serve to simplify the transformation matrix construction 
without sacrificing overall generality. 

As can be seen from Eq. (14a), the first coordinate direction in the new transformed coordinate system is 
chosen to be parallel to X 2 -X v As was mentioned earlier, the chromosome with the best fitness for the 
first objective is always placed into X v and the chromosome with the best fitness for the second objective 
is always placed into X 2 This convention is crucial for the success of the transformation algorithm, as it 
causes the first transformation coordinate to be aligned with the pareto front endpoints. Keep in mind that 
it is imperative that the first element of the P vector— the element that controls passthrough — is large 
enough so that both X, and X 2 are always passed through to the next generation without modification. It 
is also imperative for the selection algorithm to have an endpoint retention option. 

The unitary transformation matrix U is constructed column by column using 


first column 



second column 


(15b) 


(15c) 

Once G n is obtained the modification operators are applied the same as without transformation to obtain 
G /7+1 . The final G" +1 values are then obtained using an inverse transformation indicated by 

U^G^ 1 = G /?+1 (16) 

Because U is a unitary matrix, its inverse is simply its transpose, and thus, it is easily constructed. 

The gene space transformation option is controlled using the ITRAN control parameter. If ITRAN = 0 , no 
gene space transformation in implemented. If ITRAN = 1 the gene space transformation algorithm is 
activated. 


Vi = r# - u 2A u /A , Uj 2 = 

i '/ 


n m column 


/ 7— 1 


V/ ~ r i,n X U nJ u i,j > u i,n j^jj 


Geometry Parameterization 

in aerodynamic shape optimization the geometric parameterization is an important step that effectively 
connects the aerodynamic analysis routine to the optimization routine. An analytic parameterization 
suitable for wing and wing-fuselage configurations, typicai of the transonic flow regime, is now described. 
Most of the parameters have common-sense definitions in which the name itself provides the definition, 
for example, the wing leading edge radius of curvature, the wing leading edge sweep or the fuselage 
length. 

A Cartesian coordinate system ( x, y, z) is used for all configurations. The coordinate system origin is at the 
wing root leading edge for isolated wing configurations and at the fuselage nose for wing-fuselage 
configurations. A plane of symmetry ( y = 0) is inherently assumed for all configurations. The x coordinate 
is aligned with the freestream direction and is positive downstream, y is normal to the symmetry plane, 
positive to the piioi's right and z is vertical, positive upward. Aii iength parameters are nondimensionaiized 
using the wing root chord, that is, the wing chord length where it intersects the plane of symmetry. 

Isolated wino configurations — Isolated wing geometries are parameterized using N w airfoil defining 
sections or stations, each using the geometric parameterization of Sobieczky. 46 The first defining station is 
always at the wing root and the last is always at the wing tip. Each defining airfoil section is characterized by 
ten parameters. These parameters along with brief definitions are listed in Table 1 . A graphical description 
of these parameters is presented in Fig.1. All angles are measured in degrees. Once these parameters 
are specified the airfoil coordinates (z as a function of x) are determined using a polynomial of the form 

z=i^ 1/2 (17) 

n~\ 

where the coefficients a n are computed from the ten airfoil defining parameters. Two applications of Eq. 
17 are required for a complete airfoil specification — one for the upper surface coordinates and one for the 
lower surface coordinates. Six simultaneous linear equations involving rle, xup, zup, zxxup, ate + /3te/2 
and zte are solved for the upper surface coordinates, and six simultaneous linear equations involving rie, 
xio , zlo , zxxlo , ate - /3te/2 and zte’are solved for the lower surface coordinates. Because rle and zte are 
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used for both surfaces, slope continuity at the leading edge and zero thickness at the trailing edge are 
always maintained. 



AXIAL COORDINATE - x/c 


Fig. 1 Airfoil parameterization used for each wing defining station (taken from 
Sobieczky 46 ). Angles are measures in degrees. All lengths are nondimensionalized by 
wing root chord. 

Once the airfoil coordinates are constructed for each defining station, wing coordinates are constructed 
using linear lofting and then modified according to a set of planform parameters. Seven planform 
parameters are utilized at each airfoil defining station, producing a total of lxN w planform parameters. 
These parameters are also listed in Table 1 along with brief definitions. 

Certain geometric parameters have predetermined values or are not formerly used in the wing definition 
process. For example, the root chord length is always unity, and the spanwise location for the root airfoil 
defining station is always zero The tip airfoil defining station values for the dihedral angle and the leading 
edge sweep angle have no meaning, but are nevertheless included in the list for each wing 
parameterization. 


Table 1. Definitions for the airfoil section and wing planform parameters. A value for 
each parameter is required at each of the /Vowing defining stations, j=\ -,N w . 


Symbol Def in iiion * - : - r: >> W: - t: _ Symbol Definition e rr " ?■' v t 

Airfoil section parameters 

Planform parameters 

r!e y 

Leading edge radius of curvature 

c j 

Airfoil chord length 

xup 7 

x-location of upper surface max 
thickness 

d j 

Dihedral angle (deg) 

zup y 

Upper surface max thickness 

A / 

Leading edge sweep (deg) 

ZXXUp y 

Upper surface curvature at x = xup y 

'/ 

Airfoil section thickness multiplier 

XlO/ 

7 

x-location of lower surface max 
thickness 

«/ 

Wing twist (deg) 

ZlOy 

Lower surface max thickness 

*9J 

Center of twist 

ZXXlOy 

Lower surface curvature at x= xlo y 

// 

Spanwise defining station location 

ztey 

Position of trailing edge above z= 0 



ate y 

Angle between trailing edge bisector 
and z= 0 plane (deg) 



P ley 

Trailing edge included angle (deg) 




Wing-fuselage confiourations — Wings for wing-fuselage configurations are parameterized using the 
procedure as described above. Thus, only the fuselage portion of the wing-fuselage parameterization will 
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be discussed in this section. The first step is to define a base fuselage, which is analytically constructed in 

three sections. The nose section is an ellipsoid of revolution. The main body section is a right circular 

cylinder that smoothly transitions into the nose section. And the boatta.il section is a sine curve of j 

revolution smoothly connecting the main portion of the fuselage with a downstream sting. The sting — not I 

formally part of the fuselage — is a small-radius right circular cylinder that extends to the downstream f 

outflow boundary. A total of eight parameters — defined in Table 2— are required to specify this base 

fuselage. Note: The y-location of the wing root leading edge is always assumed to be on the symmetry 

plane. Thus, y R is inherently assumed to be zero and is only included for completeness. I 

Table 2. Base fuselage parameter definitions. All lengths are nondimensionalized by 

airfoil root chord. I 


Symbol 

Definition 

X H 

x-location of wing root leading edge 

Yr 

y-location of wing root leading edge 

Zr 

z-location of wing root leading edge 

X B 

Body length 

X N 

Length of ellipsoid nose 

r B 

Radius of cylindrical portion of fuselage 

X T 

Length of fuselage boattail 

r s 

Radius of fuselage sting 


Modification of the base fuselage shape is achieved using a series of Hicks-Henne bump functions. 47 A 
total of N x bump functions are distributed axially with equal spacing along the x-direction for each of N r 
circumferential stations also distributed with equal spacing in the circumferential direction. This allows a 
total of N x xN T decision variables that are used to produce incremental radius perturbations A/s to the 
base fuselage. Any position on the base fuselage, except the nose x= 0 or tail x= x B , can be changed 
using this approach. A sketch showing the definitions of these parameters, including one possible 
arrangement of A/s is shown in Fig. 2. 



Fig. 2 Sketch showing decision variable definitions for the fuselage parameterization. 
For this arrangement there are three axial stations (A/ x = 3) and four circumferential 
stations (A / r = 4). 

Masking array — Associated with each decision variable, that is, each gene, is an integer value — either one 
or zero — stored in a masking array called mask. The masking array is established at the beginning of the 
optimization computation and tells the GA which parameters to modify in the optimization process — using 
the crossover and mutation operators — and which to leave unmodified. For example, the mask values for 
Cy, y Nw , A , D Nw and y R are always zero. These parameters, as described above, are carried along with 
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each chromosome but are never used by the GA process. The masking array can also be used to limit the 
size of the design space. Only those genes with mask = 1 are utilized in the search for optimal problem 
objectives. 

Computed Results 

Area Error Norm for Two-Obiective Problems 

In order to evaluate the accuracy and efficiency of an optimization algorithm it is important to have a proper 
error norm to assess the level of convergence. For single objective problems a suitable error norm is the 
simple arithmetic difference between the current “best fitness” and the exact answer. This, of course, 
works well for model problems in which the exact answer is known. For applications in which the exact 
answer in not known a priori, it can stiii be employed as an a posteriori error computation. 

For multi-objective optimization a workable error norm is more elusive since a range of solution 
values — the so-called pareto front — is being sought. The topic of how to define easy-to-implement and 
accurate error norms for comparing one pareto front approximation to another for the same problem is 
discussed at length in Zitzler, et a!. 11 and Knowles and Corne. 18 A suitable norm encompassing the optimal 
values of each of the individual objectives is one possibility. However, such a norm would not take into 
account the trade-off regions of the pareto front. Another possibility is the attainment surface approach of 
Fonseca and Fleming. 48 In this approach a set of equally spaced sampling lines that intersect the full 
breadth of the pareto front are used. Statistical measures of goodness can then be developed based 
upon how many intersection points one pareto front has that are superior to another pareto front. 

In the present study the area between the current pareto front and the final pareto front — a numerically 
computed approximation of the finai pareto front— is used as the error norm to determine level of 
convergence. As the size of this quantity goes to zero, the current solution approaches the final solution. 

The area is computed numerically by dividing the region between the “exact” pareto front and the current 
approximation into a series of triangles as shown in Fig. 3. The area error norm is the sum of each of the 
individual triangular areas. Thus, an approximation to the exact pareto front that fails to match the final 
solution in any location will produce a non-zero value for the area error norm. 



Fig. 3 Area error norm computational process for a two-objective pareto front. 
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Isolated wing results 


The first computed result involves a two-objective optimization for an isolated wing flying in the transonic 
speed regime — AC, = 0.84 . The two objectives — iift-to-drag ratio at fixed lift and the inverse of the 
nondimensiona! structural mass — are simultaneously maximized by the previously presented MOGA 
optimization procedure. The aerodynamic Iift-to-drag ratio is evaluated for each candidate design, that is, 
for each chromosome, using the TOPS full potential solver 49 by iterating on angle of attack to obtain the 
Iift-to-drag ratio at the specified value of lift — C L = 0.45. The tolerance on the lift iteration is ±1%, but in 
most cases the error is much less. In effect, this is drag minimization at fixed lift. 

The structural mass computation utilizes an equivalent box beam model similar to that used in Oyama. 50 
Given a set of loads from the aerodynamics routine, the structures model computes the minimum-mass 
box beam that will support that load with a factor of safety of two. Shear and bending are included in this 
model but not torsion. 

The design space parameterization for isolated wings with two defining stations, N w -2 , as described 
above, consists of 34 parameters. Airfoil defining station one is at the wing root and station two is at the 
wing tip. Maximum and minimum constraint values for each of the geometric parameters, as well as 
masking array values, are given in Table. 3. Except for wing twist, all pianform geometric parameters have a 
masking array value of zero, that is, they are not modified in the optimization process. Thus, a total of 22 
genes are active within each chromosome for this optimization problem, that is, N G = 22 . 

Table 3. Maximum and minimum gene constraints along with masking array values for 
the isolated wing optimization problem. 


Parameter max ! min mask 

Parameter max min mask 

Defining airfoil section 1 

Defining airfoil section 2 

rle 1 

0.02 

0.008 

1 ' 

rle 2 

0.02 

0.008 

1 

xuPt 

0.45 

0.35 

1 

xup 2 

0.45 

0.35 

1 

zup-i 

0.07 

0.06 

1 

zup 2 

0.07 

0.06 

1 

zxxup 1 

- 0.2 

- 0.8 

1 

zxxup 2 

- 0.2 

- 0.8 

1 

xlo 1 

0.45 

0.35 

1 

xlo 2 

0.45 

0.35 

1 

zlo-i 

- 0.06 

- 0.07 

1 

zlo 2 

- 0.06 

- 0.07 

1 

zxxlo, 

0.8 

02 

1 

zxxlo 2 

0.8 

02 

1 

zte-i 

0.01 

- 0.01 

1 

zte 2 

0.01 

- 0.01 

1 

ctte! 

8.0 

0.0 

1 

ate 2 

8.0 

0.0 

1 


15.0 

6.0 

1 

ple 2 

15.0 

6.0 

1 

Pianform defining station 1 

Pianform defining station 2 

G 

1.0 

1.0 

0 

G> 

0.5 

0.5 

0 

4 

0.0 

0.0 

0 

4 1 

0.0 

0.0 

0 

At 

37.0 

37.0 

0 

CM 

< 

37.0 

37.0 

0 

4 

1.0 

1.0 

0 

h 

1.0 

1.0 

0 


3.0 

0.0 

1 

0 2 

0.0 

- 3.0 

1 

x e,i 

0.5 

0.5 

0 

X 62 

0.5 

0.5 

0 

Y\ 

0.0 

0.0 

0 

r 2 

2.4 

2.4 

0 


Results from this multi-discipline optimization problem are presented in Figs. 4-9. For all computations in 
this series of results N c = 32, /3 = 0.1, p, = 0.2, 0.2 and P=( 0.04,0.32,0.32,0.32), This P vector 

resulted in two passthrough chromosomes and ten chromosomes in each of the three remaining 
modification categories— -random average crossover, perturbation mutation and mutation. 
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Pareto front comparisons are presented in Figs. 4. Note that each comparison contains five curves that 
correspond to different ISEED initialization values used in the random number generator. There is also a 
“master” result plotted in each part of Fig. 4 that is a concatenation of all pareto front data produced for this 
case — a total of twenty curves — with only the non-dominated points retained. 




a) ISELECT = 1 , ITRAN = 0 c ) ISELECT = 3 , ITRAN = 0 



b) ISELECT = 1, ITRAN = 1 


d) iSELECT = 3, ITRAN = 1 


Fig. 4 Pareto front variation for transonic wing optimization in which the iift-to-drag 
ratio at fixed lift and the inverse of the structural mass are simultaneously maximized. 
Five solutions utilizing different seed values are plotted for each GA variation after 
500 generations, 44 = 0.84, C L = 0.45, N C = Z2 N G = 22, /3 = 0.1, 0.2, /^=0.2 and 

P= (0.04, 0.32, 0.32,0.32). 
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There are several interesting features to note in this series of pareto front comparisons: 

First of all, the variation in level of convergence across the various ISEED values is quite striking, especially 
for the cases involving greedy selection without gene space transformation, ISELECT = 1, !TRAN = 0 
(Fig. 4a). For these cases, the optimal lift-to-drag ratio point on the pareto front converges consistently in 
the number of generations allotted for this computation, regardless of the ISEED value, but the minimum 
structural mass point does not converge consistently. 

This situation is improved, somewhat, for the (SELECT = 1, ITRAN = 1 case (Fig. 4b) where the best 
structural mass endpoint values are obtained. However, this latter set of curves suffers dramatically in the 
middle portion of the pareto front where the two objectives form trade-offs with each other. 

The best results for this multi-discipline isolated-wing optimization are realized when bin selection is 
utilized, iSELECT = 3 (Figs. 4c and 4a). The pareto fronts for both sets of curves that utilize bin selection 
more closely approximate the master curve than their counterparts using greedy selection. The set of 
curves in Fig. 4c — the ISELECT = 3 , ITRAN = 0 case — are overall the best. Nevertheless, the minimum 
structural mass point for this set of computations still experiences inconsistency in convergence. 

Convergence history comparisons for each of the pareto front curves presented in Figs. 4 are presented 
in Fig. 5. Computation of the pareto front error was achieved using the area error norm described above. 
The master pareto front displayed in Figs. 4 was used as an approximation to the “exact” solution. In 
addition, for each set of convergence history curves displayed in Figs. 5 there is an average convergence 
history curve plotted (the dashed curve), which is computed by arithmetically averaging each point in the 
convergence history data file. 

Note that occasionally the convergence history error increases with generation. This usually occurs early 
in the convergence process when the pareto front is still crudely defined with a small number of discrete 
points and is explained in the series of sketches shown in Figs. 6. Figure 6a shows a hypothetical two- 
objective pareto front, both the exact curve and a crude approximation to it that might exist after the search 
has proceeded for n generations. Figures 6b and 6c show two different scenarios for what the pareto 
front might look like after n + 1 generations. 

In both scenarios a single new point has been added which neither dominates nor is dominated by any of 
the previously existing points on the pareto front. In scenario 1 the new point’s placement is less 
advantageous than that of scenario 2. Nevertheless, both new points are possible outcomes of the 
genetic algorithm search process. Despite the fact that scenario 1 represents an “enrichment” of the 
pareto front, the area error norm actually increases in going from the /7 th generation to the /?+ 1 st 
generation. This is why the area error norm occasionally jumps — either up or down — during solution 
evolution. 

There is a moderately large amount of variance in each of the convergence history curves displayed in Fig. 
5. This emphasizes the point that genetic algorithms are stochastically based. Comparison of results, 
especially convergence efficiency results, must be made with the proper amount of statistical averaging in 
place. Results from each MOGA variation demonstrate this conclusion. 
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a) ISELECT = 1, ITRAN = 0 


c) ISELECT = 3, ITRAN = 0 




b) ISELECT = 1, ITRAN = 1 d) ISELECT = 3, ITRAN = 1 


Fig. 5 Convergence efficiency comparisons for transonic wing optimization in which 
the lift-to-drag ratio at fixed lift and the inverse of the structural mass are 
simultaneously maximized. Five solutions utilizing different seed values are plotted for 
each GA variation, M x = 0.84, £^=0.45, A/ c = 32 N a - 22, /3 = 0.1, /^=0.2, /%=0.2 and 
P- (0.04, 0.32,0.32,0.32). 


The averaged convergence history efficiencies for each of the four MOGA variations studied herein are 
presented on the same set of axes in Fig. 7. Note that the two variations that utilize bin selection are both 
faster than their greedy selection counterparts, achieving more than an order of magnitude reduction in 
the area error norm after 500 generations. The gene-space transformation variation produced mixed 
results, being slightly faster when used with greedy selection and slightly slower when used with bin 
selection. 
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OBJECTIVE 1 OBJECTIVE 1 OBJECTIVE 1 • j 

a) Baseline, pareto front b) Pareto front after n + 1 c) Pareto Troni arier n + i ; 

after n generations. generations (scenario i). generations (scenario 2). 

Fig. 6 Sketch of pareto front convergence history process showing two different 
advancement scenarios. | 



0 100 200 300 400 500 

GENERATION NUMBER 

Fig. 7 Average convergence history curves for the four GA variations piotted in Figs. 4 

and 5, M„ =0.84, C L = 0.45, N c = 32 N a = 22, p = 0.1, /3, = 0.2, /*, = 0.2 and 

P= (0.04, 0.32, 0.32, 0.32). 

Computed results for the two-objective isoiated-wing optimization problem taken from the endpoints of a 
suitable ISELECT = 3, ITRAN = 0 pareto front are displayed in Figs. 8-10. Figure 8a shows the upper 
surface Mach number contours for the minimum drag point (at fixed lift), and Fig. 8b shows the same view 
for the minimum structural mass point. Note the extreme suction pressure (high Mach number) upstream 
of the shock wave in Fig. 8b. This is caused by a large amount of flow expansion on the forward part of the 
wing and results in a much stronger shock wave than the solution depicted in Fig. 8a. 

This situation can be viewed in a more quantitative way by looking at cross-sectional plots of the pressure 
coefficient, which are displayed in Figs. 9 for three spanwise stations — y!b = 0.21, 0.51, 0.78. In each plot 
the pressure coefficient for the two pareto front endpoints — minimum drag and minimum structural 
mass— -are displayed. From this series of figures it can be seen that the minimum drag case (solid curves) 
has a greatly reduced shock strength relative to the minimum structural mass case (dashed curves). 
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Both solutions have the same total lift. Despite this, the spanwise distribution of lift is different between 
the two cases. This can be seen from Figs. 9 by carefully examining the areas circumscribed by each 

mum structure! mass solution {css hod curvos) contains 
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more area inboard and iess outboard relative to the minimum drag point. This tends to move the center of 
pressure inboard for the minimum structural mass solution and thus, reduces the moment arm through 
which the wing lift acts. This, in turn, reduces the wing root bending moment and allows for a lighter 
structure. 



Fig. 8 Upper surface Mach number contours for a transonic wing optimization in which 
the drag at fixed lift and the structural mass are simultaneously minimized, ISELECT = 3, 
ITRAN = 0, = 0.84, b7 £ =0.45, /V~ = 32 N a = 22, /3 = 0.1, /z = 0.2, /x, = 0.2 and 

P- (0.04, 0.32. 0.32, 0.32). 
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Fig. 9 Pressure coefficient distributions at selected spanwise stations for the wing 
optimization problem presented in Fig. 8 showing differences between the two pareto 
front endpoints, iSELECT = 3, iTRAN = 0, M x = 0. 84, C =0.45, N r = 32 /V- = 22, p = 0.1 ? 
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Another important attribute in this multi-discipline optimization problem is wing thickness— especially at 
the wing root. Selected airfoil cross-sectional profiles are presented in Fig. 10 at ylb = 0.0, 0.48, 0.98 
(wing root, mid-span, wing tip). Profiles for both pareto front endpoints — minimum drag (solid curves) and 
minimum structural mass (dashed curves)— are presented. Note that the vertical axis had been greatly 
expanded to facilitate viewing of the results, and that each profile is displayed in rigged position with wing 
twist included. For this set of computations the minimum-drag root-airfoil thickness is the minimum allowed 
by the geometric constraints and the minimum-structural-mass root-airfoil thickness is the maximum 
allowed. This is as expected. A reduced wing-root thickness helps reduce drag by reducing shock 
strength and thus reduces wave drag. An increased wing-root thickness helps to reduce structural mass 
by more efficiently supporting the wing root bending moment with a box beam that has increased depth. 



NORMALIZED AXIAL DISTANCE - x/c 

Fig. 10 Selected airfoil profiles in rigged positions — at y/b=0.0, 0.48 and 0.98 — for the 
transonic wing optimization problem of Fig. 8 showing differences between the two 
pareto front endpoints, (note the expanded z/c scale), ISELECT = 3, ITRAN = 0, 

4/. = 0.84, C L = 0.45 , A/ c = 32 Afe = 22, £ = 0.1, tf=0.2, #,=0.2 and #=(0.04,0.32,0.32,0.32). 

Wino-fuselaoe results 

The next set of computed results involves a two-objective optimization for a wing-fuselage configuration 
flying in the transonic flow regime — =0.84. The two objectives — lift-to-drag ratio at fixed lift and the 
configuration's volume — are simultaneously maximized by the previously presented MOGA optimization 
procedure. As before, the aerodynamic lift-to-drag ratio is evaluated for each candidate design, that is, for 
each chromosome, using the TOPS full potential solver 49 by iterating on angle of attack to obtain the lift-to- 
drag ratio at the specified value of lift — C L = 0.45 . 

The wing-fuselage volume computation utilizes a straightforward algorithm that divides the fuselage into a 
series of cross-sections, computes the area of each cross-section using a simple triangularization and 
then computes the volume by integrating along the fuselage using trapezoid integration. The wing 
volume is computed using a similar algorithm that is deactivated for those wing cross-sections that lie 
inside the fuselage. 

There are several reasons lift-to-drag ratio (at fixed lift) and wing-fuselage volume have been chosen as 
the two objective functions for the present optimization problem. First of all there is a clear-cut conflict 
between these two objectives. When the volume is maximized the lift-to-drag ratio suffers and vice versa. 
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Thus, we have a classical engineering trade-off problem. How the MOGA performs on such a problem with 
real-world objectives is the key aspect of this study, not the development of new aeronautical engineering 

r'Anoonfc 

^ 1 lOOpiO . 

In addition, the vehicle volume was chosen as an objective because it is an easy quantity to understand. 
The direct computation of maximum volume is easy to perform, thus providing a simple check for the 
results obtained by the GA — at least for one endpoint on the pareto front. This problem is interesting in 
that it is complicated by transonic-flow shock-wave-induced non linearities and by a large number of 
decision variables that have large sensitivity variations. 

The present wing-fuselage parameterization with two wing defining stations, N w = 2, as described above, 
consists of 66 parameters. Airfoil defining station one is at the wing root and station two is at the wing tip. 
There are 24 A r values used to define the fuselage — six equally-spaced axial stations, each with four 
equally-spaced meridinal locations. The axial positions are located at a 7 £=0.71,1.43, 2.1 4, 2.86, 3.57, 4.29 
where the fuselage length is fixed at x B = 5.Q. The meridinal positions are located at 
(, o = -90°, - 30°, 30°, 90° where <p = -90° corresponds to the fuselage keel line and q> = 90° corresponds 
to the crown line. 

Maximum and minimum constraint values for each of the geometric parameters, as well as masking array 
values, are given in Table. 4. As can be seen, 43 of the 66 genes that make up each chromosome have 
masking array values of one, and 23 have values of zero, that is, only 43 genes are actually modified 
during the optimization process— N G = 43 . Note also that the max-min values associated with each of the 
A r values are not symmetric about zero. Along the keel and crown lines, for example, any A r value 
averaged between the max and min constraints, will be positive, and along the other meridinal stations 
(along the fuselage side), the same averaging will produce negative values. The max and min A r 
constraints were chosen in this way io keep the vertical extent of the fuselage from becoming too small. A 
difficulty arises in the wing-fuselage line of intersection computation if the intersection line contacts the 
symmetry plane. This results in fuselage designs that have depth-to-width ratios that generally exceed 
unity — sometimes by a substantial amount. 

The first wing-fuselage optimization results are presented in Figs. 11-14. For these results the freestream 
Mach number is 0.84 and the lift coefficient Is held fixed at 0.45. The tolerance on the lift iteration is ±1%, 
but in most cases the error is much less. The second binning selection algorithm with endpoint retention 
is used, ISELECT = 4, and the gene-space transformation procedure is turned off, ITRAN = 0. Both of 
the mutation probabilities are fixed at 0.2 and P= (0.04, 0.32, 0.32, 0.32) . 

This optimization uses 34 chromosomes for each generation. The first element of the P vector is set such 
that two passthrough chromosomes are chosen during each generation — always the pareto front 
endpoints — -the current minimum drag and the maximum volume chromosomes. Thus, 32 modified 
chromosomes require new function evaluations during each generation. These function evaluation 
computations are efficiently performed simultaneously on 32 processors of a parallel computer. The other 
three elements of the P vector are set to evenly divide the remaining chromosomes between the three 
remaining modification operators — random average crossover, perturbation mutation and mutation (at 
least as close to even as possible). 

Development of the pareto front as a function of generation number is displayed in Fig. 11. Note that the 
maximum volume point on the pareto front, which occurs at a value of 1 .42, converges quickly, achieving 
96% of the theoretical maximum in as few as 50 generations, while the minimum drag point is slower in 
convergence. This characteristic, displayed in Fig. 11 for one ISEED value, occurred for every ISEED 
value utilized and is due to several factors. 
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Table 4. Maximum and minimum gene constraints along with masking array values for 
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First of all, the volume objective depends on most of the genes in a simple way. For example, if one of the 
A r genes is increased, while all other genes are held fixed, the volume increases in a proportionate 
amount throughout the range of that gene’s variation. This behavior exists no matter what value the other 
genes possess. In addition, the maximum volume values of many genes in the present wing-fuselage 
parameterization — certainly all of the fuselage A r values — are at their upper constraint values. Finding an 
optimum on the design space boundary seems to be easier than finding an optimum within the design 
space interior. 

The drag objective depends on the various genes in a more complex manner than that of vehicle volume 
— especially for this wing-fuselage problem. Interrelated and simultaneous changes in several genes are 
typically required in certain regions of the design space to achieve improvement in the drag objective. For 
example, when z R is changed the position of the wing-fuselage juncture changes. To achieve optimal 
drag performance the wing-fuselage area ruling needs to be simultaneously adjusted, and this requires 
coordinated changes in various elements of the Ar array. Which elements have to be changed and by 
how much depends on the value of z R . 
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The slow convergence of the drag pareto front endpoint for the present wing-fuseiage problem, relative 
to the second objective, was not observed in the previous problem. In fact, the second objective for the 
previous isolated-wing optimization problem — minimization of the structural mass — was slower than that of 
the drag. There are two reasons for this apparent inconsistency. First of all, the complex interdependence 
of the drag objective on wing placement and the ensuing fuselage area ruling issues do not arise in the 
isolated-wing problem. Secondly, the dependence of the structural mass .on the gene value distribution 
for the isolated-wing problem is every bit as complicated as for the drag objective — perhaps more so. 
That’s because the structural mass depends on both the wing thickness distribution, especially near the 
wing root, as well as the wing pressure distribution — the latter quantity setting the wing’s center of 
pressure and thus the moment arm for the wing root bending moment. 



Fig. 11 Pareto front convergence and solution structure for a wing-body two-point 
optimization involving drag minimization and volume maximization both at constant lift, 
C L = 0.45, =0.84,”lSELECT = 4, ITRAN = 0, N a -2, 7^ = 34, N a = 43, 0 = 0.1, p, =0.2, 

P 2 ~ 0.2 and P = (0.04, 0.32,0.32,0.32). Wing sectional pressure distributions are presented 
at three positions along the pareto front, A) maximum volume, B) intermediate, C) 
minimum drag. 

Also displayed in Fig. 11 are several airfoil pressure distributions for three different solutions 
corresponding to three different locations along the pareto front, A) maximum volume, B) intermediate 
and C) minimum drag. For each solution two airfoil pressure distributions are presented — the first inboard 
near the wing-fuselage juncture (y! b-= 0.2) and the second near mid-semi-span {ylb- 0.6). The most 
obvious feature from these results is the change in shock strength as the pareto front is traversed. The 
maximum volume solution (A) contains strong shocks, both on the upper and lower surfaces of the wing, 
while the minimum drag solution (C) contains mild shocks on only the wing’s upper surface. 

Upper surface Mach number contours for each of the three solutions presented in Fig. 1 1 are displayed in 
Figs. 12. This series of figures shows the shock wave pattern variation on the upper surface of the wing- 
fuselage configuration as the solution transitions from the maximum volume pareto front endpoint to the 
minimum drag endpoint. For the maximum voiume endpoint there are strong shock waves, not only on the 
wing, as seen in Fig. 11, but also on the fuselage near the nose and tail. The fuselage shock waves are 
caused by the GA’s selection of minimum values for the x N and x r genes. These values maximize the 
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vehicle volume, but cause shock waves to form at both the fuselage nose and tail due to excessive over- 
expansion. 



A) Max volume solution 



5) Intermediate solution 



C) Min drag solution 


Fig. 12 Surface Mach number contours for the wlng/bcdy optimization problem 
presented in Fig. 11 at three positions along the pareto front. 
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noticeable initial change is associated with the fuselage nose and tail fairings. By choosing larger values 
for x N and x T , the drag is decreased without sacrificing a iarge amount of volume. This can be seen by 
examining the intermediate solution in Fig. 12 and noting its position on the pareto front in Fig. 1 1 . Fortner 
reductions in drag are achieved via the complex task of fuselage area ruling in the vicinity of the wing- 
fuselage juncture. Such a result is observed by looking at the minimum drag solution in Fig. 12. More on 
this, including a comparison of selected fuselage cross-sections will be presented subsequently. 

Pressure distributions along the fuselage for the three solutions described in Fig. 11 are presented in 
Figs. 13. Resuits are plotted along three fuselage meridina! stations — the keel line (<? = - 90 c ), the side, 
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MOGA optimizer has chosen a high mounted wing, regardless of the solution’s position along the pareto 
front. More will be presented on this aspect subsequently. Because of this the side fuselage pressure 
distribution (cp = 0°) is plotted along a continuous meridina! line that lies below the wing for each of these 
three solutions. In addition, Mach number contours as viewed from the fuselage’s side are also displayed 
in Figs. 13. 


As can be seen from Fig. 13a — and previously discussed in conjunction with Hg. 12 — the expansions, 
especially at the fuselage nose and tail, are quite large and lead to strong shock waves. As the fairings at 
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fhe shock wave induced on the fuselage by wing-fuselage interference can also be seen in Fig. 
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c) Minimum drag solution. 

Fig. 13 Fuselage pressure distributions for the wing/bcdy optimization problem 
presented in Fig. 11 at three positions aiong the pareto front. 


By comparing each part of Figs. 13 with the corresponding part in Figs. 12. it can be seen that the depth cf 
the fuselage, as measured from keel to crown is much larger than twice the half-width for any given 
longitudinal station. This is a direct result of the choices made in setting up the constraints for the A r array, 
as was already discussed in the fuselage parameterization section., and produces fuselage cross-sections 
that are elongated in the z-direction. This situation Is presented in a more quantitative way in Fig. 14. 



A series of fuselage cross-sections are presented in Fig. 14 for the two-objective optimization problem 
displayed in Fig. 11. For each cross-section there are two curves plotted. The solid curve corresponds to 
the solution obtained at the minimum-drag pareto front endpoint, and the dashed curve corresponds to 
the solution obiained at the maximum-volume pareto front endpoint. For this solution the wing root 
leading edge is placed at {x Rt y R ,z R ) = (1.5, 0.0, 0.1 5). Note that ** = 1.5 and y R = 0.0 are values set at 
the beginning of the computation and are held fixed, and that z R = 0.1 5 is a value selected by the MOGA 
optimization for this particular chromosome from the range, -0.1 5 <2* <0.1 5 (see Table 4). The only 
fuselage cross-section in Fig. 14 that actually intersects with the wing is the one at */c=1.54. As seen 
from Fig. 1 4, the fuselage cross-section is drawn before the wing intersection is included. 

From this set of plots it is easy to see how the fuselage was area ruled in the vicinity of the wing-fuselage 
juncture to obtain the final amount of drag reduction Volume is subtracted along the upper portion of the 
fuselage near where the wing intersects, creating “pear-shaped” cross-sections in this fuselage location. 
Every cross-section from the maximum volume solution contains more area than the corresponding 
minimum drag cross-section. Each of the A r values for the maximum volume solution is individually 
maximized. 



Fig. 14 Selected fuselage cross-sections for the wing/fuseiage optimization problem 
of Fig. 11 showing differences between the two pareto front endpoints, £^ = 0.45, 
At =0.84, ISELECT = 4, ITRAN = 0, N 0 = 2, /V c = 34, Afc = 43, >3 = 0.1, F\ = 0.2, Afe=0.2 and 
A 1 = (0.04, 0.32, 0.32, 0.32). 

Convergence efficiency for the two-objective wing-fuselage optimization problem just discussed — drag 
minimization and volume maximization (both at fixed lift) — is presented in Figs. 15 and 16. Figure 15 
shows the previously discussed area error norm plotted versus the number of function evaluations for 
three different convergence histories involving three different population sizes — 34, 66 and 130 
chromosomes, respectively. Each curve is averaged over two ISEED values in an attempt to eliminate 
some of the statistical variation. As can be seen each convergence history is nearly the same, indicating 
that convergence is not a function of population size. This is compatible with the results presented in 
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Holst and Pulliam 15 where a number of model problems were used to evaluate GA convergence for two- 
objective optimization. 

For this set of computations the number of processors used to perform the function evaluations is equal 
to the number of chromosomes requiring function evaluations. There are always two passthrough 
chromosomes. The remaining chromosomes — those needing modification — are arbitrarily divided into 
three groups of equal size — or as close to equal as possible. One group is modified using crossover, the 
second using perturbation mutation and the third using global mutation. Since, two passthrough 
chromosomes are utilized for each population, the number of processors required to perform function 
evaluations is 32, 64 and 128, respectively. Using an increased number of chromosomes to define a 
given population allows an increased number of processors to be used for the optimization problem and 
thus, improves turn-around efficiency. 


Figure 1 6 presents a comparison of the six pareto fronts computed from the convergence efficiency study 
presented in Fig. 15. These results are statistically similar to previous results presented in this section and 
suggest that there is no significant effect of population size on the final solution. Care should be taken in 
evaluating this conclusion, as the number of solutions needed to construct a proper statistical average has 
not been utilized in this case due to computer time limitations. 



Fig. 15 Pareto front convergence history 
comparisons for three population sizes, 
each averaged over two (SEED values, 
wing/fuselage optimization problem. 

C L = 0.45, AC =0.84, ISELECT = 4, 

ITR AN = 0 , N 0 = 2, N G = 4 3, 0 = 0.1, #=0.2, 
£> = 0.2. 



LIFT TO DRAG RATIO 


Fig. 16 Pareto front comparisons for 
three population sizes, each with two 
different ISEED values, wing/fuselage 
optimization problem, £^=0.45, AC =0.84, 
ISELECT = 4, ITRAN = 0 , N a = 2, /V G = 43, 
0 = 0 . 1 , # = 0 . 2 , #= 0 . 2 . 


It is of interest to study how the individual genes vary along the pareto front. Selected genes from the 
wing-fuselage optimization problem described in Fig. 11 are plotted along the pareto front in Figs. 17 and 
18. In particular, they are plotted as a function of the pareto front’s first objective, lift-to-drag ratio. Each 
gene is normalized and then plotted so that it varies from 0 to 1 . For many of the curves the discrete points 
are plotted along with a smoothed result, in some cases, where the number of curves on a single set of 
axes is large, only the smoothed curve is plotted. 
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a) .r-location of airfoil maximum 
thickness on the lower surface at the 
wing tip, GENE^ = xlo 2 . 



b) Airfoil trailing edge angle at wing root 
and tip, GENE 7 = ate! and GENE 14 =ate 2 . 


c) Airfoil twist at wing root and tip, 

GENE 15 =0, and GENE 16 =0 2 . 



d) Vertical wing position ( GENE 17 = z R ), 
fuselage nose length ( GENE 18 = x N ) and 
fuselage boattail length (GENE ig = x r ). 


Fig. 17 Normalized gene distributions along the pareto front from selected positions 
within each chromosome, wing/fusel age optimization problem. ^=0.45, At =0.84, 
ISELECT = 4, ITRAN = 0 , N 0 = Z, /V c = 34, N G = 43, jS = 0.1, A =0.2, /%=0.2, 

P- (0.02, 0.33, 0.33, 0.32). 

The first general observation is that a significant amount of statistical variation exists in many of the genes 
as they vary across the pareto front. Some genes wildly oscillate without any discernible pattern. The x- 
location of the wing-tip lower-surface maximum thickness — xlo 2 , gene number 1 1 — is such an example. Of 
the 43 available genes, xlo 2 contains the most variation. This gene is plotted in Fig. 17a and is seen to 
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favor the maximum value in its allowable range. Nevertheless, it takes on many other values in a seemingly 
haphazard fashion. 

The reason for this behavior is straightforward. This gene has little impact on either objective being 
optimized and thus has little impact on. the pareto front. Any value is generally acceptable at any location 
along the pareto front. To test this observation several chromosomes were selected from the pareto front 
for which xlo 2 was against one of its constraints. The value of xlo 2 within that chromosome was then moved 
to the opposite constraint without changing any other gene value, and a new solution was computed. The 
resulting design-space point moved little and was still a member of the pareto front. Ideally, such a 
parameter shouldn’t be utilized as a gene in the optimization process. In reality it is difficult to know a priori 
which genes will behave in this way and which will not. 

Other parameters piotted in Fig. 17 include the trailing edge angle at the wing root and tip (Fig. t7b), the 
airfoil twist at the wing root and tip (Fig. 17c), and the vertical wing mounting position, the nose length and 
the boattail length (Fig. 17d). All of these gene parameters have reasonably controlled variations across 
the pareto front. 

Figure 18 presents the variation of each of the A r parameters for the wing-fuselage optimization 
problem — all 24 of them — across the pareto front. As already discussed, the A r parameters are organized 
into four meridinal stations, each with six longitudinally distributed values. Each meridinal station is plotted 
on a separate set of axes in Fig. 18. Figures 18a, 18b, 18c and 18d show the A r longitudinal variation 
along <p = -90°, - 30°, 30° and 90° , respectively. 


As expected, all A r values at the maximum volume endpoint on the pareto front lie at their maximum 
constraint limits. Moving away from this endpoint, many of the A r values decrease from their maximum 
constraint values and attain intermediate values. This is especially true for the A r values along the 0 = 30° 
and 0 = 90° meridinal stations, where reductions in the various A r values have been determined by the 
MOGA optimizer to achieve optimal fuselage area ruling in the vicinity of the wing-fuselage juncture. 

It is of interest to study the effect of wing placement on the overall wing-fuselage optimization problem. In 
the cases presented above, the vertical wing position, z R , is constrained to lie between -0.15 and 0.15, 
while x R and y R are fixed at 1.5 and 0.0, respectively. As a result of the optimization, the value of z R is 
chosen to lie at or near its upper constraint value, z R = 0. 15, over the entire pareto front (see the variation 
of gene 17 plotted along the pareto front in Fig. 17d). In other words, a high wing intersection is favored 
by the optimization regardless of position on the pareto front — at ieast in the context of the present 
design space construction. 


Figures 1 9 and 20 present results for a wing-fuselage optimization in which the wing intersection is 
constrained to be low mounted, that is, the value of z R is constrained to be -0.15. All other design space 
attributes and constraints are the same as in the previous wing-fuselage computations. This is 
implemented by setting the maximum constraint for z R equal to the minimum constraint, z R = -0.1 5, and 
by setting the corresponding z R masking array value to zero. Thus, for this probiem, N G - 42, instead of 
the previous value of 43. Note also that N c = 66 for this series of computations, that is, 64 processors 
were used in the optimization. 


ITS 3. Pi c! tWO 


Figure 1 S displays two computed pareto fronts for the original wing-fuselage optimization prt 
pareto fronts for the optimization that constrains the wing to be low-mounted. The two curves plotted for 
each case correspond to two different ISEED values. As can be seen, the low mounted wing results are 
less efficient in the vicinity of the minimum-drag pareto front endpoint. The low-mounted case achieves 
essentially the same level of drag, but only when vehicle volume is reduced. Note also that this effect is 
large in comparison to the statistical variation that exists between the two ISEED computations from each 
solution category. The effect of z R on the maximum volume pareto front endpoint is negligible. 
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Fig. 18 Normalized gene distributions along the pareto front from the fuselage portion 
of the parameterization, C L =0.45, M. = 0.84, ISELECT = 4, ITRAN = 0, N 0 = 2, Afc = 34, 
N g = 43, /3 = 0.1, =0.2, P2= 0.2 , P= (0.02, 0.33, 0.33, 0.32). 

Also displayed in Fig. 19 are two circular symbols showing discrete chromosomes from each different 
category of pareto front that most cioseiy match a Lift-to-drag ratio of 13. Selected fuselage cross sections 
taken from these two points are compared in Fig. 20. A totai of five longitudinal stations are presented with 
x/c= 0.49, 1 .54, 2.59, 3.64, 4.69. The solid curve corresponds to the original case in which the wing 
vertical mounting position was allowed to float. (Although for this chromosome the value of z R that the 
MOGA selected was 0.15— a high mounted wing). The dashed curve corresponds to the case in which 
the wing was constrained to be low mounted. Despite being at the same L7D location on the pareto front, 
the two area rulings are quite different — quite different at each longitudinal station. Close examination 




shows how the MOGA modified the area ruling to account for the shifted wing mounting position, taking 
volume from the upper fuselage for the high-mounted case and from the lower fuselage for the low 
mounted case. It is also obvious that the fuselage associated with the low-mounted wing contains a 
smaller amount of volume as is consistent with Fig. 19. 
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Fig. 19 Pareto front comparisons for two different wing-fuselage optimization cases: a) 
baseline in which wing position is allowed to float ( -0.15 < z R < 0.15) and b) wing fixed in 
the low-wing position (z R = -0.15). Each computation performed with two ISEED values, 

C L = 0.45 , 44 =0.84, (SELECT = 4, ITRAN = 0 , N a = 2, /V c = 66, Afe = 43/42, 0 = 0.1, #=0.2, 
#>=0.2, F= (0.02, 0.33, 0.33, 0.32). 
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Fig. 20 Selected fuselage cross-sections for the two wing/fuselage optimization 
problems of Fig. 19 showing differences between the two solutions at LiD = 13: a) 
baseline in which wing position is allowed to float (- 0.1 5 <z R < 0.1 5 ) and b) wing fixed in 
the low-wing position (^ = -0.15), C L = 0.45, 44=0.84, ISELECT = 4, ITRAN = 0, /V a = 2, 
N c = 66, N g = 43/ 42, 0 = 0.1, # = 0.2, #>=0.2, #=(0.02,0.33,0.33,0.32). 
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Conclusions 


A multi-objective genetic algorithm (MOGA) optimization procedure, with several selection algorithm 
options and a new gene-space transformation procedure, has been presented. It uses real-number 
encoding to represent all design space decision variables as genes and populations of fixed size to go 
from generation to generation. Four modification operators are utilized to advance from one generation to 
the next. They include passthrough, random average crossover, perturbation mutation and mutation. The 
standard output for this approach is a pareto front, which includes the best solutions from each objective, 
as well as a range of optimal “tradeoff” solutions in between. 

The MOGA optimization procedure was utilized to determine a number of multi-objective optimal solutions 
for two problems. The first involved a transonic wing lift-to-drag maximization (drag minimization) coupled 
with the simultaneous minimization of structural mass (ail at fixed lift). The second involved a transonic 
wing-fuselage lift-to-drag maximization (drag minimization) coupled with the simultaneous maximization of 
vehicle volume (all at fixed lift). In every case the GA optimization process was convergent. However, some 
of the GA variations studied converged more rapidly than others. 

Of the two selection algorithms compared — greedy selection and bin selection — the latter produced the 
most efficient convergence across all of the problems studied. The gene-space transformation procedure 
produced mixed results in the area of GA convergence efficiency. In some cases it was moderately faster, 
for other cases moderately slower. 

The new MOGA utilizes a masking array, which allows any gene (or set of genes) to be eliminated from the 
optimization process as a modifiable decision variable. This allows determination of the effect of a single 
gene (or gene subset) on the resultant pareto front. 

A limited parallelization efficiency study was performed involving the use of 32 to 128 processors. 
Population sizes ranging from 34 to 130 chromosomes were used. The present MOGA is essentially 
embarrassingly parallel with regard to the number of chromosomes used for each generation, and the 
number of chromosomes did not affect convergence efficiency of the MOGA scheme. Thus, the 
turnaround time dramatically decreased with increasing number of processors. 
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